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Abstract — We present a hierarchical model predictive control 
approach for large-scale systems based on dual decomposition. 
The proposed scheme allows coupling in both dynamics and 
constraints between the subsystems and generates a primal 
feasible solution within a finite number of iterations, using 
primal averaging and a constraint tightening approach. The 
primal update is performed in a distributed way and does 
not require exact solutions, while the dual problem uses an 
approximate subgradient method. Stability of the scheme is 
established using bounded suboptimality. 

I. Introduction 

Coordination and control of interacting subsystems is an 
essential requirement for optimal operation and enforcement 
of critical operational constraints in large-scale industrial 
processes and infrastructure systems [1]. Model Predictive 
Control (MPC) has become the method of choice when de- 
signing control systems for such applications [2]-[4], due to 
its ability to handle important process constraints explicitly. 
MPC relies on solving finite-time optimal control problems 
repeatedly online, which may become prohibitive for large- 
scale systems due to the problem size or communication 
constraints. Recent efforts have been focusing on how to 
decompose the underlying optimization problem in order to 
arrive at a distributed or hierarchical control system that 
can be implemented under the prescribed computational and 
communication limitations [5], [6]. One common way to 
decompose an MPC problem with coupled dynamics or con- 
straints is to use dual decomposition methods [7]-[9], which 
typically lead to iterative algorithms (in either a distributed or 
a hierarchical framework) that converge to feasible solutions 
only asymptotically. Implementing such approaches within 
each MPC update period can be problematic for some 
applications. 

Recently, we have presented a dual decomposition scheme 
for solving large-scale MPC problems with coupling in both 
dynamics and constraints, where primal feasible solutions 
can be obtained even after a finite number of iterations [10]. 
In the current paper we present a novel method that is 
motivated by the use of constraint tightening in robust 
MPC [11], along with a primal averaging scheme and dis- 
tributed Jacobi optimization. Since an exact optimum of the 

The authors are with Delft Center for Systems and Control, 
Delft University of Technology, Delft, The Netherlands {m.d.doan, 
t. keviczky, b. deschutter}@tudelf t . nl 



Lagrangian is not assumed to be computable in finitely many 
iterations, an approximate scheme is needed for solving the 
MPC optimization problem at each time step. We present a 
solution approach that requires a nested two-layer iteration 
structure and the sharing of a few crucial parameters in 
a hierarchical fashion. The proposed framework guarantees 
primal feasible solutions and MPC stability using a finite 
number of iterations with bounded suboptimality. 

The paper is organized as follows. In Section |IIJ we 
describe the MPC optimization problem and its tightened 
version, which will be used to guarantee feasibility of the 
original problem even with a suboptimal primal solution. 
Section [III] describes the main elements of the algorithm 
used to solve the dual version of the tightened optimiza- 
tion problem: the approximate subgradient method and the 
distributed Jacobi updates. In Section IIVI we show that 
the primal average solution generated by the approximate 
subgradient algorithm is a feasible solution of the original 
optimization problem, and that the cost function decreases 
through the MPC updates. This allows it to be used as a 
Lyapunov function for showing closed-loop MPC stability. 
Section [VT] concludes the paper and outlines future research. 

II. Problem description 

A. MPC problem 

We consider M interconnected subsystems with coupled 
discrete-time linear time-invariant dynamics: 

M 

•<•:.., Y^-^-'i • i = i,...,M (i) 

and the corresponding centralized state-space model: 

x k+ i = Ax k + Bu k (2) 

with Xk = [{xi) T (xir...(xf) T r, Uk = 

[{ul) T {ul) T ...(uf) T ] T , A = -AikjL { L... ; M} and 
B = [Sy]ij e {i,...,M}- 

The MPC problem at time step t is formed using a convex 
cost function and convex constraints: 

t+N-l , v 

min ^ ixl"Qx k + ulRu k \ + xJ +N Px t +N (3) 
k=t ^ ' 



s.t. 4+1= ^ A«4 + s««i, 

i = 1, . . . , M, k = t, . . . , t + N - 1 (4) 

x k £X,k = t+l,...,t + N-l (5) 

x t+N e x f c X (6) 

u k €U,k = t,...,t + N -1 (7) 

ai6fii,i = l M, fc = i, ...,t + N-l (8) 

a* = z(t) G A" (9) 

where u = [uf, . . . , uf+ iv _ 1 ] T , x = [xf +1 , . . . , a;f +JV ] T , 
the matrices Q, P, and i? are block-diagonal and positive 
definite, the constraint sets U, X and Xf are polytopes and 
have nonempty interiors, and each local constraint set fij is 
a hyperbox. Each subsystem i is assigned a neighborhood, 
denoted J\f % , containing subsystems that have direct dynami- 
cal interactions with subsystem i, including itself. The initial 
state xt is the current state at time step t. 

As U, X and Xf are polytopes, the constraints © and 
© are represented by linear inequalities. Moreover, the state 
vector x is affinely dependent on u. Hence, we can eliminate 
state variables Xt+i, • • ■ , %t+N and transform the constraints 
©, ©, and (O into linear inequalities of the input variable 
u. Eliminating the state variables in ©-©I leads to an 
optimization problem in the following form: 

/; = min f(u,x t ) (10) 

U 

s.t. g(u,x t )<0 (11) 
ueO (12) 

where / and g = [g\, , . . ,g m ] T are convex functions, and 
fi = YitLi ^2 with each fi; = Jl^o 1 ^« * s a hyperbox. 
Note that f(u,x t ) > 0,Vu ^ 0,x t ^ 0, due to the positive 
definiteness of Q, P, and R. 

We will use (ut,x t ) to denote a feasible solution gener- 
ated by the controller for problem ©-© at time step t. 
This solution is required to be feasible but not necessarily 
optimal. We will make use of the following assumptions: 

Assumption 2.1: There exists a block-diagonal feedback 
gain K such that the matrix A + BK is Schur (i.e., a 
decentralized stabilizing control law for the unconstrained 
aggregate system). 

Assumption 2.2: The terminal constraint set Xf is posi- 
tively invariant for the closed-loop Xk+i = (A + BK)xk 
(x e wk(Xf) =>■ (A + BK)x G int(Af f )). 

Assumption 2.3: The Slater condition holds for problem 
(fT0]i-(fT2l. i.e., there exists a vector that satisfies strict in- 
equality constraints [12]. It is also assumed that prior to each 
time step t, a Slater vector u t is available, such that 

gj(vL t ,x t ) < Q,j = l,...,m (13) 
Remark 2.4: Since g(u,xt) < has a nonempty interior, 
so do its components gj(u,xt) < 0,j — l,...,m. Hence, 
there will always be a vector that satisfies the Slater condition 
(fT3l . In fact, we will only need to find the Slater vector Qo 
for the first time step, which can be computed off-line. In 
Section IV- Al we will show that a new Slater vector can then 
be obtained for each t > 1, using Assumption 12.21 



Assumption 2.5: At each time step t, the following holds 

f{u t -\,x t -i) - f(u u x t ) > xJ^Qxt-i +uf_ l Ru t -i 

(14) 

For later reference, we define A f > which can be 
computed before time step t as follows: 

A f = xf^Qxt-i + uJ^Rut-x (15) 
Remark 2.6: Assumption 12.51 is often satisfied with an 

appropriate terminal penalty matrix P. A method to construct 

a block-diagonal P with a given decentralized stabilizing 

control law is provided in [13]. 

Assumption 2. 7: For each xt G X, the Euclidean norm of 

g(u,x t ) is bounded: 

Lt>||ff(u,a;t)||2,Vuen (16) 
Remark 2.8: In the first time step, with given xq, we 
can find Lq by evaluating ||g(u, xq)||2 at the vertices of O, 
the maximum will then satisfy (fTol l for t = 0, due to the 
convexity of g and f2. For the subsequent time steps, we 
will present a simple method to update L t in Section IV-BI 

B. The tightened problem 

We will not solve problem (TTOt — (fT2b directly. Instead, we 
will make use of an iterative algorithm based on a tightened 
version of (TTOt— (fT2b . Consider the tightened constraint: 

g'(n,xt) =g(u,x t ) + l m c t < (17) 

with g'(u,x t ) = [g[,. . . ,g' m } T , < c t < 

mirij = i m {—gj(u t , Xt)}, and l m the column vector 

with every entry equal to 1, Due to ( fT3l . we have 

max {g'Ju t ,Xt)} = max {gj(u t ,x t )} + Ct < (18) 

j—l,....rn J j— l,...,m 

Hence g'j(u t ,xt) < 0,j = 1, . . . , m. Moreover, using (TToT l 
and the triangle inequality of the 2-norm, we will get L' t — 
L t +c t as the norm bound for g', i.e. L' t > \\g'(vL, Xt)\\2, Vu G 
O. Note that L' t implicitly depends on x t , as u t and ct are 
updated based on the current state xt- 

Using the tightened constraint (I17t . we formulate the 
tightened problem: 

/r=min f(u,x t ) (19) 
u 

s.t. g'{u,x t )<0 (20) 

u g n (2i) 

Remark 2.9: Only the coupled constraints (flTT l are tight- 
ened, while the local input constraints (fT2l are unchanged. 
The Slater condition also holds for the tightened problem 
(fT9T>— (l2TT>. with Of being the Slater vector. 

III. The proposed optimization algorithm 

Our objective is to calculate a feasible solution for problem 
©-(O using a method that is favorable for distributed 
computation. The main idea is to use dual decomposition for 
the tightened problem (fT9b — d2TT > instead of the original one, 
such that after a finite number of iterations the constraint 
violations in the tightened problem will be less than the 
difference between the tightened and the original constraints. 



Thus, even after a finite number of iterations, we will obtain 
a primal feasible solution for the original MPC optimization 
problem. 

A. The dual problem 

We will tackle the dual problem of ([T9}-(|2"T1>. in order to 
deal with coupled constraint g'(u,Xt) < in a distributed 
way. In this section, we define the dual problem and its 
subgradient. For simplicity, in this section the dependence 
of functions on the initial condition xt is not indicated 
explicitly. 

The Lagrangian of problem (fT9b — (f2TT> is defined as: 

£'(u,/i) = /(u) + /iV(u) (22) 

in which u 6 fi, \i e R+. 

The dual function for (QjO-lED: 

4(p) =min£ / (u,/z) (23) 

is a concave function on R™, and it is non-smooth when / 
and g 1 are not strictly convex functions [12]. 

Given the assumption that Slater condition holds for ([T9l>- 
(l2Tl i. duality theory [12] shows that: 

?r = /r (24) 

with (ft = maxjjgiT. <z'(/i) and / t '* the minimum of (fT9b — 

(ED. 

Thanks to this result, instead of minimizing the primal 
problem, we may maximize the dual problem, which is often 
more amenable to decomposition due to simpler constraints. 
Since we may not have the gradient of q' in all points of 
R™, we will use a method based on the subgradient. 

Definition 3.1: A vector d is called a subgradient of a 
convex function / over X at the point x £ X if: 

f(v) > f(x) + (y - xfd, VyeX (25) 

The set of all subgradients of / at the point x is called 
the subdifferential of / at x, denoted df(x). 

For each Lagrange multiplier p G R™, first assume we 
have u(/Z) = argmin ue ^ C(u, p). Then a subgradient of 
the dual function is directly available, since [12]: 

< + (M - £)V (ufa)), e R™ (26) 

In case an optimum of the Lagrangian is not attained due 
to termination of the optimization algorithm after a finite 
number of steps, a value u(p) that satisfies 

£'(u(/2),/2) < mm£'(u,/Z) +6 (27) 
will lead to the following inequality: 

<Z'( M ) < </(/2) + 6 + (n - p) T g'(u(p)), Vai e R" 1 (28) 

where </(u(/i)) is called <5-subgradient of the dual function 
<7 at the point /2. The set of all (S-subgradients of q at p is 
called (5-subdifferential of q at /2. 

This means we do not have to look for a subgradient (or 
<5-subgradient) of the dual function, it is available by just 
evaluating the constraint function at the primal value u(p) 
(or u (/!)). 



B. The main algorithm 

We organize our algorithm for solving (TTOb— (fT2b at time 
step t in a nested iteration of an outer and inner loop. The 
main procedure is described as follows: 

Algorithm 3.2: Approximate subgradient method with 
nested Jacobi iterations 

1) Given a Slater vector u t of (fT0b-(fT2l. determine c t and 
construct the tightened problem (fT9b— (f2TT). 

2) Determine step size at and suboptimality e t , see later 
in Section InFcT] 

3) Determine k t (the sufficient number of outer itera- 
tions), see later in Section IIII-C.2I 

4) Outer loop: Set p,W = q . i m . For k = 0, ... , k, find 
u (fc) jjU (fc+i) such that . 

£'(u (fe) ,/i (fe) ) <xmn£'(u,pW)+e t (29) 

M (fc+1) =^{^ (fc) +«*d (fe) } (30) 

where Vvl™ denotes the projection onto the nonnegative 
orthant, d^ = g'(u^\x t ). 
Inner loop: 

• Determine pk (the sufficient number of inner iter- 
ations), see later in Section Ull-D. 11 

• Solve problem ( 1291 in a distributed way with 
a Jacobi algorithm. For p = 0, ...,pk, every 
subsystem i computes: 

11*0+1) =arg min C'^ip), . . . ,u i - 1 (p),u i ! 

u l+1 (p),...,u M (p), M W) (31) 

where fli is the local constraint set for control 
variables of subsystem i. 
. Define u« = [u 1 {p k ) T , . . . , u M {p k ) T } T , which 
is guaranteed to satisfy ( |29b . 

5) Compute u (fet) =^Ez=o u(0 ' take u « = ^ as the 
solution of (fT0b-(fl2l. 

Remark 3.3: Algorithm [372] is suitable for implementation 
in a hierarchical fashion where the main computations occur 
in the Jacobi iterations and are executed in parallel by local 
controllers, while the updates of dual variables and common 
parameters are carried out by a higher-level coordinating 
controller. In the inner loop, each subsystem only needs to 
communicate with its neighbors, which will be discussed in 
Section IIV-AI This algorithm is also amenable to implemen- 
tation in distributed settings, where there are communication 
links available to help determine and propagate the common 
parameters k t , and p k . 

In the following sections, we will describe in detail 
how the computations are derived, and what the resulting 
properties are. 

C. Outer loop: Approximate subgradient method 

The outer loop at iteration k uses an approximate sub- 
gradient method. The primal average sequence ir^ = 
£ 2~2i=o U ^ nas f°ll owm g properties: 



?'(fl 



For k > 1 



< 



1 



7t 



-[/(u t ,x t ) -qf] 



<*tL^ 
27* 



(32) 



e* (33) 



where g' + denotes the constraint violation, i.e. 
max{</,0 • l m }. The proof of (l32l can be found in 
and the proof of d33l is given in Appendix IVII-AI 

7) Determining at and St- Using the lower bound of the 
cost reduction (TBI and the upper bound of the suboptimality 
(l33l for the tightened problem ([T9ll-(|2"T1>. we will choose at 
and s t such that f(u t ,x t ) < /(u t _i, x t _i). 

The step size at and suboptimality e t should satisfy: 



.9 



/+ 



[14], 



a t L' t 



■ e* < At 



(34) 



where At is defined in (TT3T >. and LJ is the norm bound for 
cf . This condition allows us to show the decreasing property 
of the cost function in problem ©-(O, which can then be 
used as a Lyapunov function. 

Note that a larger at will lead to a smaller number of outer 
iterations, while a larger et will lead to a smaller number of 
inner iterations. For the remainder of the paper we choose 
their values according to 



At 
2 



(35) 



(36) 



2) Determining k t : Using the constraint violation bound 
\, we will choose k t such that at the end of the algorithm, 
we will get a feasible solution for problem (IToTi— (IT2Ti. which 
is the average of primal iterates generated by 



The subgradient iteration 
1, . . . , kt, with the integer 



1 



(=0 



(37) 



is performed for k = 



1 



a t Ct \lt 



2j t 



+ a t L' 



(38) 



defined a priori, where ["•] is the ceiling operator which gives 
the closest integer equal to or above a real value, 7 t = 



=l,„.,m{-g'i(V-t,Xt)} 



and Ut is the Slater vector of (fT9]>— (|2TT> . 



i{-9i(ut,x t )} 



D. Inner loop: Jacobi method 

The inner iteration (f3TT > performs parallel local optimiza- 
tions based on a standard Jacobi distributed optimization 
method for a convex function C'(u,fi^) over a Cartesian 
product, as described in [15, Section 3.3]. In order to find 



the sufficient stopping condition of this Jacobi iteration, we 
need to characterize the convergence rate of this algorithm. In 
the following, we summarize the condition for convergence 
of the Jacobi iteration, noting that C'(vi,fi^) is a strongly 
convex quadratic function with respect to u. 

Proposition 3.4: Suppose the following condition holds: 



(39) 



where with i,j S {1, . . . , M} denotes a submatrix of 
the Hessian H of £' w.r.t. u, containing entries of H in 
rows belonging to subsystem i and columns belonging to 
subsystem j, A m i n means the smalleast eigenvalue, and a 
denotes the maximum singular value. 

Then 3</> £ (0, 1) such that the aggregate solution of the 
Jacobi iteration OTb satisfies: 

\Hp) -u*|| a < M0 p max||u 4 (O) -u"|| 2 , Vp > 1 (40) 

i 

where u* = argmin ue f2 £'(u, /i ( - fc - ) ), and u ! * is the compo- 
nent of subsystem i in u*. 

We provide a proof for Proposition l3.4l in Appendix lVII-BI 
Remark 3.5: This proposition provides a linear conver- 
gence rate of the Jacobi iteration, under the condition of 
weak dynamical couplings between subsystems. For the sake 
of illustrating condition (l39l . let all subsystems have the 
same number of inputs. Consequently, Hij is a square and 
symmetric matrix for each pair hence the maximum 

singular value a (Hij) equals to the maximum eigenvalue. 
Inequality (l39l thus reads: 

\nin{Hu) > ^ A max (LTjj), Vi 

which implies that the couplings represented by H are small 
in comparison with each local cost. 

Remark 3.6: Note that the strong convexity of £' and the 
condition (l39l are required only for the convergence rate 
result of the Jacobi iteration in which C is a quadratic 
function. Extensions to other types of systems, where the 
Lagrangian can be solved with bounded suboptimality, are 
immediate. In such cases we simply need to replace the 
Jacobi iteration with the new algorithm in the inner loop, 
while the outer loop will remain intact. 

1) Determining p^: As £'(u, •) is continuously differen- 
tiable in a closed bounded set ft, it is Lipschitz continuous. 

Suppose we know the Lipschitz constant A of £'(u, ■) over 
SI, i.e. for any u 1 , u 2 G Jl the following inequality holds: 



|£V.M (fc) )-£'(uV fc) )l|2 < AHu 1 



u 2 || 2 



(41) 



Taking u 1 = u(pk) and u 2 = u* in (PiTT i. and combining 
it with (l40l . we obtain: 



\C'(u(p k ),^)-mmC'(u,^)\\ 2 < A\\u(p k ) - u*\\, 



< AM(j) pk max||u 4 (0) - u 4 



(42) 



For each i e {1, . . . , M}, let D. L denote the diameter of the 
set fij w.r.t. the Euclidean norm, so we have ||u*(0) — U l * ||a < 
Di. Hence the relation (l42l can be further simplified as 

£'(u(p fe ), M (fe) )<min£'(u,/iW) 



uer2 



h.M(fr k max A (43) 



Based on (1431 . in order to use u(pk) as the solution 
that satisfies (|29T i, we choose the smallest integer such 
that AMf' max,; A < e t : 



Pk 



log 



AM max, D, 



(44) 



IV. Properties of the algorithm 



A. Distributed Jacobi algorithm with guaranteed conver- 
gence 

The computations in the inner loop can be executed by 
subsystems in parallel. Let us define an r-step extended 
neighborhood of a subsystem i, denoted by Af£, as the set 
containing all subsystems that can influence subsystem i 
within r successive time steps. AC is the union of subsystem 
indices in the neighborhoods of all subsystems in N^_i'. 



AC 



u & 



(45) 



where Af{ = TV 4 . We can see that in order to get update in- 
formation in the Jacobi iterations, each subsystem i needs to 
communicate only with subsystems in A/]y-_ 1 , where N is the 
prediction horizon. This set includes all other subsystems that 
couple with i in the problem (TTOb — (fT~2l > after eliminating the 
state variables. This communication requirement indicates 
that we will benefit from communication reduction when the 
number of subsystems M is much larger than the horizon N, 
and the coupling structure is sparse. 

Assume that the weak coupling condition ([39l l holds, then 
after pf. iterations as computed by d44l i. the Jacobi algorithm 
generates a solution u(" = n{pk) that satisfies (|29l l in the 
outer loop. 

B. Feasible primal solution 

Proposition 4.1: Suppose Assumptions 12.11 and 12.31 hold. 
Construct g 1 as in ( fTTI i. a t as in (l35l >. Let the outer loop 
(g5)-(ED) with fj,W = ■ l m be iterated for k = 0, . . . , k t . 
Then u < - fet - ) is a feasible solution of ([ToTi-([T2"li. where u^*-* is 
the primal average, computed by (|37T >. 

Proof: With a finite number of k t iterations (l32t reads as 



< ■= 

2 k t Ci t 



7t 
2lt 



It 



a t L' t 



(46) 



Moreover, the dual function q' t is a concave function, there- 
fore q'* > q'(0,x t ). Recall that f{u,x t ) > 0,Vu ^ 0,x t ^ 
0, thus q'(0,x t ) = min uer! f(u,x t ) + • lf n g'{u,x t ) = 
min ug0 f(n,x t ) > 0, thus 



x t 



< r~~ {—Si^t.xt) 



+ -7— + a t L t 
tit 



(47) 



Combining d47b with ( [38] ), and noticing that k t and c t are 
all positive lead to 



^g'(u Ckt \x t ) <ct, j = l 



(48) 

/ \ " " " - ; / ' / ■ ../ ....... w . (49) 

^^(fi^,^) <0, j = l,...,m (50) 

where the last inequality implies that u' fe *^ is a fea- 
sible solution of problem (fT0]i-(fT2l. due to ct < 



-l,...,m{-g 3 (Ut,X t )}. 



C. Closed-loop stability 



□ 



Proposition 4.2: Suppose Assumptions 12.31 12.51 and 12.7 
hold. Then the solution u^ fc *' ) generated by Algorithm 13.2 



satisfies the following inequality: 

f{u u x t ) < f(u t -i,x t -i), VteZ+ (51) 
Proof: Using (l33l l and (l34l >. and noting that fi^ = 0, we 
obtain: 



V / 2k t a t 2 

(52) 



Notice that u t is also a feasible solution of (fT~9T> — (f2TT > 
(due to the way we construct the tightened problem: u t 
still belongs to the interior of the tightened constraint set), 
while //* is the optimal cost value of this problem. As a 
consequence, 



r; <f(ut,x t ) 



(53) 



Combining (|52]l, d53), and (O, and noting that u t = u (fc ^ 
leads to: 



f(vL t ,x t ) < /(u t _i,x t _i), Vi e Z 4 



(54) 



□ 

Note that besides the decreasing property of f(ut,Xt), all 
the other conditions for Lyapunov stability of MPC [16] 
are satisfied. Therefore, Proposition 14.21 leads to closed- 
loop MPC stability, where the cost function f(u tl xt) is a 
Lyapunov candidate function. 

V. Realization of the assumptions 

In this section, we discuss the method to update the Slater 
vector and the constraint norm bound for each time step, 
implying that Assumptions 12.31 and 12 .7 1 are only necessary in 
the first time step (t = 0). 



A. Updating the Slater vector 

Lemma 5.1: Suppose Assumption ^. 2l holds. Let u t be the 
solution of the MPC problem d3j-(|9]l at time step t, computed 
by Algorithm [372] Then iit+i constructed by shifting Ut one 
step ahead and adding ut+N = Kxt+N, is a Slater vector 
for constraint ([TO at time step t + 1. 

Proof: Note that based on Proposition 14.11 is a 

feasible solution of problem dTOb — (fT2T>. Moreover, the strict 
inequality ( T50T > means that u^*- 1 is in the interior of the 
constraint set of Oil-®. This also yields: 

x t+ N S mt(Xf) (55) 

Moreover, due to Assumption 12.21 we have (^4 + 
BK)x t +N G int(Af). This means that if we use ut+N = 
Kx t +N, then the next state is also in the interior of the 
terminal constraint set Xf. Note that hi and X do not change 
when problem (O-© is shifted from t to t + 1, hence all 
the inputs of Ut+x and their subsequent states are in the 
interior of the corresponding constraint sets. Therefore, Ut+x 
as constructed at step 5 of Algorithm 13.21 is a Slater vector 
for the constraint (flTT i at time step t + 1. □ 

This means we can use Ut+x = u*+i as the qualifying 
Slater vector for Assumption 12.31 at time step t + 1. 

B. Updating the constraint norm bound 

In our general problem setup, g(u, x) is composed of affine 
functions over u and x, and thus can be written compactly 

as 

g(u, x) = Ex + 9u + r (56) 

with constant matrices 2, O and vector r. Then for each 
Xt-i, x t , and u e fi, the following holds: 

g(u, x t ) = g(u, x t -x) + S(x t - x t -i) 

ll.9(u,xt)|| 2 < \\g(u,x t _x)h + l|S(a: t -a; t _i)||a (57) 

In order to find a bound L t for g(u,x t ) in each t > 1 
step, we assume to have the constraint norm bound available 
from the previous step: 

L t -i > ||g(u,x t „i)|| 2 ,Vu G n (58) 

Hence, combining the above inequalities a norm bound 
update for g(u,x t ) can be obtained as: 

L t = L t _i + ||H(a:t-i t _i)||2 (59) 

VI. Conclusions 

We have presented a constraint tightening approach for 
solving an MPC optimization problem with guaranteed feasi- 
bility and stability after a finite number of iterations. The new 
method is applicable to large-scale systems with coupling 
in dynamics and constraints, and the solution is based on 
approximate subgradient and Jacobi iterative methods, which 
facilitate implementation in a hierarchical or distributed way. 
Future extensions of this scheme include a posteriori choice 
of the solution by comparing the cost functions associated 
with the Slater vector u t and the primal average u ( - fct ^ in a 
distributed way. 
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VII. Appendix 
A. Proof of the upper bound on the cost function (133b 

This proof is an extension of the proof of Proposition 3(b) 
in [14], the main difference being the incorporation of the 
suboptimality e t in the update of the primal variable d29l >. 

Using the convexity of the cost function, we have: 

\ i=o / 1=0 



\ J2 (/(n«) + (m W )V(« W )) - I J> (0 ) VO^) 



1=0 



1=0 



(60) 



Note that £'(nW,/iW) = f /(u«) + <?'(uW)V° J and 

£'(u«,M (,) ) < min£'(u«, M (i) ) +£* = ?'(/i (0 ) 

VZ < k 
(61) 

Combining the two inequalities above, we then have: 

/(u W ) < I ]T + * - J E(m w )V(u«) 



<?r+^-^E(^ ) ) TdW 



(62) 



1=0 



where dS 1 ' — g'(u^), and the last inequality is due to 

q' t * > 

Using the expression of squared sum: 

\\^ +1) \\l<y {l) +a t d^\\l 

= ||M W ||2 + 2at(A* W ) r d (,) + 11^111 (63) 

we have: 



-(^^^^(ll^llI-llM^lli + afll^ll 



for I = 0, . . . , k - 1. 

Summing side by side for / = 0, . . . , k — 1, we get: 



(64) 
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(65) 



Linking ( 1621 and ( I65l l. we then have 



/(u«)<gf +e t 
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|^ (0) || 2 - ||/i (fe) || 2 
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i=0 



< 



'7/ 



|| M (°) || 2 , a t Lr 



2fca t 



(66) 



in which we get the last inequality by using L' t as the norm 
bound for all g'(u®), 1 = 0, k - 1. 

Finally, with the Slater condition, there is no primal-dual 
gap, i.e. q' t * = / t * (cf. d24l>). hence: 



□ 



5. Proof of the convergence result of the Jacobi iteration 
( Proposition I3.4D 

According to Proposition 3.10 in [15, Chapter 3], the 
Jacobi algorithm has a linear convergence w.r.t. the block- 
maximum norm, as defined below: 

Definition 7.1: For each vector x = [xj , . . . , xj^] with 
Xi £ M™% given a norm || ■ ||j for each i, the block-maximum 
norm based on II ■ ||, is defined as: 



x b -m = max ^ 



Definition 7.2: With any matrix A G 



(67) 
, we asso- 



ciate the induced matrix norm of the block-maximum norm: 
\\AWn = maxM^ = max \\Ax\U (68) 

x=i0 \\X\jj |M|j-=l 

In this paper, we use the Euclidean norm as the default 
basis for block-maximum norm, i.e. || • ||j = || • 1 1 2 , Vi . 

Proposition 3.10 in [15, Chapter 3] states that u(p) gen- 
erated by OTI ) will converge to the optimizer of C'(u,x t ) 
with linear convergence rate w.r.t. block-maximum norm 
(i.e. \\u( P ) - u*|| b _ m < 0*||u(O) - u*|| b _ m , with u* = 
argmin u £'(u, x t ) and <f> G [0, 1)) if there exists a positive 
scalar 7 such that the mapping R : l~i n> K" u , defined by 
i?(u) = u — 7V u £'(u, Xt), is a contraction w.r.t. the block- 
maximum norm. 

Our focus now is to derive the condition such that R(u) 
is a contraction mapping. 

Note that since /(u, Xt) is a quadratic function, and 
g'(u,x t ) contains only linear functions, the function 
C'(u, Xt) is also a quadratic function w.r.t. u, hence it can 
be written as: 



£'(u, x t ) = u T Hu + b T u + c 



(69) 



where H is a symmetric, positive definite matrix, b is a 
constant vector and c is a constant scalar. 

In order to derive the condition for R(u) to be a contrac- 
tion mapping, we will make use of Proposition 1.10 in [15, 
Chapter 3], stating that: 

If / : R n » 1 — ^ R n " is continuously differentiable and there 
exists a scalar <\> G [0,1) such that 

||7-7Gr 1 (V i F i (u)) r || ii + E|| 7 Gr 1 (V j F i (u)) T || ii <^ 

Vu G 

(70) 



defined with each com- 

-r-l 



then the mapping T : fl n- 

ponent i G {1,...,M} by T;(u) = Uj - jG^ 1 F(a) is a 
contraction with respect to the block-maximum norm. 

The mapping T(u) will become the mapping i?(u) if we 
take d = ,Vi and F(u) = V u £'(u,x t ) = 2H» + b. 
With such choice, and evaluating the induced matrix norm 
d68l l in (f70l . the condition for contraction mapping of R(u) 
is to find <p G [0, 1) such that: 



i» 2 



•53|27#tf|| 2 <&Vi 



(71) 



where Hj with i,j G {1,...,M} denotes the submatrix 
of H, containing entries at rows belonging to subsystem i 
and columns belonging to subsystem j. Note that the matrix 
inside the first induced matrix norm is a square, symmetric 
matrix, while the matrices ify are generally not symmetric, 
depending on the number of variables of each subsystem. 
The scalar <\> £ [0, 1) is also the modulus of the contraction. 

Using the properties of eigenvalue and singular value of 
matrices, we transform (F7Tb into the following inequality: 

max|27A(iJ lt ) - l\+2 1 ^a{H l0 ) < <j>,Vi (72) 

where A means eigenvalue, and a denotes the maximum 
singular value. 

In order to find 7 > and <fr G [0, 1) satisfying d72l . we 
need: 



Note that the closer of tf> to 0, the faster the aggregate 
update u(p) converges to the optimizer of the Lagrange 
function. 

In order to get the convergence rate w.r.t. the Euclidean 
norm, we will need to link from the Euclidean norm to the 
block-maximum norm: 

M 

IW2 < V ||x l || 2 < Mmax\\x% = M|[a;|| b . m (80) 

* — ' i 

i=l 

Hence, the convergence rate of Jacobi iteration OTb w.r.t. 
the Euclidean norm is: 

||»(p) -u*||a < M0 p max||u 4 (O) -u"|| 2 , Vp > 1 (81) 

i 

□ 



max|27A(iJ Jl ) - 1| + 27^ a (Hij) < l,Vi (73) 



2 7 A 

max 



,Vi (74) 



1 - 2j\ min (H u ) + 2 7 E^ a(Hj) < 1 

«J 7 < V (W^«) + Ej* j v , (75) 

[ A ro in (Ha) > J2j& a ( H ij) 

The first inequality of (l75T l shows how to choose 7, while 
the second inequality of ( l75l l needs to be satisfied by the 
problem structure, which implies there are weak dynamical 
couplings between subsystems. 

In summary, the mapping R(u) satisfies (l70l and thus is 
a contraction mapping if the following conditions hold: 

1) For all i: 



Amin (H u )>J2°(Hij) 

2) The coefficient 7 is chosen such that: 

1 

7 < 



(76) 



(77) 



So, when condition d76l ) is satisfied and with 7 chosen by 
( fTTl i. we can define G (0, 1) as: 

cj) = max t max j 27(A max (fljj) + ^(Hj)) - 1, 

l-27(A min (^)-^a(^))| 1 



(78) 



This is the modulus of the contraction R(u), and also 
acts as the coefficient of the linear convergence rate of the 
Jacobi iteration d3Tl >, which means: 



||u(p)-u*|| b . m <<^||u(0)-u*||b. m , Vp>l (79) 
where u* = argmin uG f2 C'(a,Xt). 



